rm(list = ls())

# Note that this script takes 3.5 hours to run because it does
# 37,830 randomization inference tests with 1,000 permutations each
# skip to line 78 to use the pre-stored permutation results.

wd <- ".../Replication/"
setwd(wd)

# Load/install packages --
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  foreign, 
  ggplot2, 
  estimatr,
  texreg,
  xtable,
  fastDummies,
  sandwich,
  dplyr,
  janitor, 
  gridExtra,
  gsheet,
  zoo,
  interflex,
  lubridate,
  tidyverse,
  stringi,
  readxl,
  ri2,
  modelsummary,
  ggpubr
)

options(scipen=999)

# Note: given the random component inherent to the randomization tests performed
# by ri(), you may observe small differences in the p-values of the randomization
# tests performed using that function.

# load data
lottery <- read.csv("Data/baseline database.csv")
data <- subset(lottery, lottery$group!="Liberado por eleccion")

# Appendix F --------------------------------------------------------------------
#Results Iteratively Excluding Observations

# exclude two lottery winners at a time
data$id <- rownames(data)

permutations <- expand.grid(data$id,data$id)
permutations <- subset(permutations, permutations$Var1!=permutations$Var2)


# Note: the iteration of randomization inference over permutations takes about 3.5 hours
# Skip to line 77 for pre-stored results 

coef_loo2 = c()
p_loo2 = c()

for(i in 1:nrow(permutations)){
  
  d.temp = subset(data, data$id!=permutations$Var1[i] & data$id!=permutations$Var2[i])
  declare.temp <- declare_ra(N = nrow(d.temp), m=sum(d.temp$lottery_winner)) 
  
  ri_excluding <- summary(conduct_ri(formula = imprisoned_dummy ~ lottery_winner,
                                     declaration = declare.temp, sharp_hypothesis = 0, 
                                     assignment = "lottery_winner",
                                     data = d.temp, sims = 1000))
  
  coef_loo2[i] = ri_excluding[2]
  p_loo2[i] = ri_excluding[3]
  
  print(i)
  
}

# Skip the iteration and use the pre-stored results instead 
load("Data/Figure A8_simulated results.Rda")

min(unlist(coef_loo2))
max(unlist(p_loo2))


pdf("Output/FigureA8.pdf",height=5, width=10)
par(mfrow=c(1,2))
hist(unlist(coef_loo2),
     main = "Leave-two-out coefficients",
     xlab = "Estimated coefficient",
     xlim = c(0,.15),
     breaks = 10)
hist(unlist(p_loo2),
     main = "Leave-two-out exact p-values",
     xlab = "Exact p-value",
     xlim = c(0,.15),
     breaks = 100)
dev.off()
